Polarization Charge Distribution in Gapped Graphene 
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We study the distribution of vacuum polarization charge induced by a Coulomb impurity in 
massive graphene. By analytically computing the polarization function, we show that the charge 
density is distributed in space in a non-trivial fashion, and on a characteristic length-scale set by 
the effective Compton wavelength. The density crosses over from a logarithmic behavior below this 
scale, to a power law variation above it. Our results in the continuum limit are confirmed by explicit 
diagonalization of the corresponding tight-binding model on a finite-size lattice. Electron-electron 
interaction effects are also discussed. 

PACS numbers: 81.05.Uw, 73.43.Cd 



INTRODUCTION AND FORMULATION OF 
THE PROBLEM 



Over the course of the past year the behavior of 
graphene in the presence of a strong external Coulomb 
field was analyzed in considerable detail. 1 < 2 < 3 > 4 < 5 < 7 This 
problem is important, notably, for our understand- 
ing of electronic transport in the presence of charged 
impurities. 8 In addition, since the effective coupling con- 
stant in graphene, a = e 2 /(Hv), can be rather large (in 
vacuum a = 2.2), the exploration of features uniquely 
associated with the strong-field regime Za ~ 1 (where 
Z stands for the strength of the external Coulomb 
field: V(r) = Ze/r) becomes a realistic prospect. In 
this context, it was found that above a critical value 
(Za) c = 1/2 characteristic resonances appear in the en- 
ergy spectrum 1 ' 2 and the vacuum polarization density 
varies as p(r) ~ 1/r 2 . In the subcritical regime (Za < 
1/2), on the other hand, the polarization charge is con- 
centrated around the Coulomb center in such a way that 
one obtains p(r) oc S(r), within the continuum approx- 
imation for the electron dynamics. The physics around 
the critical point (Za) c is a peculiar masslcss realization 
of the more "conventional" vacuum charging behavior in 
massive quantum electrodynamics (QED), 6 which also 
seems possible in graphene under certain conditions. 7 

The drive to explore unconventional behavior at strong 
coupling (Za ~ 1) has been fueled, so far, by theoreti- 
cal progress only. Experiments in which K + ions are 
deposited in a controlled way onto graphene show the 
behavior one expects from the theory of scattering of 
Dirac fermions by a Coulomb field, 9 ' 10 but only in the 
undercritical regime, perhaps as expected for such low 
value of Z. It is also clear that under current experi- 
mental conditions a < 1 due to dielectric screening by 
the substrate, and additional screening is provided by 
the presence of water layers around the samples. 11 This 
conspires to significantly increase the effective dielectric 
constant with a concomitant decrease of a, thus mak- 
ing the subcritical regime Za < (Za) c the relevant one 
for present-day experiments, and also accounting for the 
feeble signatures of interaction effects. It is therefore nat- 
ural to ask how the characteristic features of the under- 



critical regime manifest themselves in the vacuum polar- 
ization and screening properties. In the strict massless 
limit the polarization charge density is concentrated at 
the potential source, p(r) cx ZaS(r), and non-trivial spa- 
tial variation can only occur due to additional interaction 
effects. 3 These are expected to be small, and we will see 
that, perturbatively, they read p(r) ~ Za 2 /r 2 . 

In the present work we explore another source of den- 
sity variation caused by the presence of a finite "mass" 
m or, equivalcntly, an energy gap A = 2mv 2 in the elec- 
tronic spectrum. There are at least three sources of a 
gap in graphene. Firstly, it has been realized recently 12 
that epitaxially grown graphene exhibits a substantial 
gap (A = 0.26 eV) due to the breaking of the sublattice 
symmetry by the substrate. Graphene suspended above 
a graphite substrate also has a small gap A « 10 meV. 13 
Secondly, spin-orbit coupling leads to a gap, albeit of 
much smaller magnitude: A so ~ 10 -3 meV. 14 Finally, 
real mesoscopic samples can have an effective gap gener- 
ated by their finite size, which scales as A ~ 1/L. We find 
that the polarization density of massive Dirac fermions 
displays characteristic behavior, controlled by the effec- 
tive "Compton" wavelength Ac = h/(mv). Our main 
results are summarized in Figs. 3 and 4: for distances 
r < Ac, the density variation is logarithmic, crossing 
over to a power law tail at r > Ac. It should be pos- 
sible to explore this behavior with modern experimental 
techniques for detection of local density variation. 13 ' 15,16 

Our starting point is a 2D system of massive Dirac 
electrons in the presence of a Coulomb impurity with 
effective charge Ze. The Hamiltonian is: 



Ze 2 

H = —ih ver ■ V + mv 0-3 

r 



(1) 



where v is the Fermi velocity and Ci are the Pauli ma- 
trices. The induced vacuum polarization charge is calcu- 
lated in linear response according to 



/>(q)=ZeV(q)II(q > 0), V(q) 



2ita 



(2) 



where a = e 2 /(hv) and II(q, 0) is the static polarization 
function. This equation is schematically represented by 
the diagram in Fig. 1. 
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FIG. 1: Graphical representation of Eq. (2). 



Unless specified otherwise, we use units in which v = 
h = 1, the electron charge is — e (e > 0) and, for conve- 
nience, we measure all charges (p, Q) in units of e (the 
sign is meaningful though). The chemical potential, //, 
is assumed to be in the gap, < m (further discussion 
appears later). 

In the rest of the paper, we first compute the polar- 
ization operator for massive graphcnc (Section II), which 
we use in Section III for the calculation of the induced 
vacuum polarization charge in a weak Coulomb field. In 
Section IV we compare the obtained behavior with re- 
sults based on exact diagonalization studies on a finite- 
size lattice. Section V discusses the influence of weak 
electron-election interactions on the polarization charge, 
and Section VI contains our conclusions. 



II. POLARIZATION FUNCTION FOR MASSIVE 
DIRAC PARTICLES 

The polarization function is computed in the standard 
way as 

U(q,u>) = -iNj2 J ^Tr{G(k,u)G(k + d,u + oj)}, 

k 

(3) 

where the trace is over the Pauli matrices and V = 4 ac- 
counts for the valley and spin degeneracies. The Green's 
function at finite mass is given by 



G(k,i/) = 



m (T3 



"I 



(4) 



Using a more symmetric 3- vector notation, (q, qg 
(q,w), (k,ko) = (k,v), fc 2 = k 2 -^ 2 



/cq, etc., we obtain 



n(q,go) 



d 3 k fc (fc + qa) + k ■ (k + q) + 



(2tt) 3 (fc 2 +m 2 ) 



[(k + qf 



(5) 

Technically, an exact analytical evaluation of (5) be- 
comes possible if one treats the frequency and momen- 
tum integrations on equal footing, as would happen in a 
Lorentz invariant situation, and was done for the mass- 
less case. 17 However in this way one encounters a linearly 
divergent piece, since it is clear that at large momenta 
Eq. (5) leads to J A d 3 k/k 2 ~ A, where A is the covari- 
ant ultraviolet cut-off. This procedure therefore requires 
"regularization" , i.e. subtraction of the infinite contri- 
bution. The regularization procedure, however, yields 



the correct result because in a non-relativistic situation, 
when the energy (kg) integration is performed first in the 
interval (— oo, oo), the resulting momentum integration is 
ultraviolet convergent (and cut-off independent). 

Restoring the original notation, we obtain the final ex- 
act expression for the polarization 



n(q,w) = - 



where 




1 - 



4m 



arctan 



2m/ 



(6) 



(7) 



Unlike QED, the polarization function of graphene is not 
covariant. We have confirmed this result by direct nu- 
merical integration of (3). A comparison of the two in 
the static case (oj = 0) is shown in Fig. 2, where the corre- 
spondence between the numerical calculation and Eq. (6) 
is exact. 

At finite frequency, dynamical properties such as the 
longitudinal conductivity of gapped graphcnc can be de- 
rived directly from Eq. (6). For the real part of the 
conductivity a m {ui) one can use the standard formula, 



18 



for 



0. We read- 



a m (uj) = -e 2 {oj/q 2 ) ImII(q, 
ily extract the imaginary part of the dynamical polariza- 
tion, 



iqr I 4m 2 
Imn(q, lo) = I 1 



4g 



u-Vhl 2 + 4m 2 ) , 



where we define 
This leads to 

CO 



q = yjoj 2 - |q| 2 . 



= 1 



4m' 



(oj — 2m) 



(8) 



(9) 



(10) 



where <7q is the conductivity of masslcss graphene, pre- 
dicted to be cr = e 2 /(4/i). 8 Eq. (10) implies that, at 
the edge oj = 2m, the conductivity increases by a fac- 
tor of two, a m = 2(7q, and decreases for higher frequen- 
cies, approaching a m — <7q (oj ^> m). For a gap of 
A = 2m ?s 260 meV, 12 we estimate that the characteris- 
tic frequencies where the enhancement should be observ- 
able are u ~ A ~ 2 x 10 3 cm -1 , which are typical fre- 
quencies in infrared spectroscopy. 19 Similar effects were 
previously discussed in magnetotransport (in particular 
when a gap is opened in strong magnetic field). 20 



III. INDUCED VACUUM POLARIZATION 
CHARGE 

In the following discussion, we assume oj = (q = |q|), 
and denote II(q) = II(q, 0). First we extract the behavior 
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FIG. 2: (Color online) Plot of the function F(m/q), defined 
as: II(q,ct; = 0) = —qF(m/q), where q = |q|. 



in some limiting cases. In the limit of large momenta we 
have 



whereas in the opposite limit 

,;2 



n(q) 



37rm 



m 



m 



< 1 



(11) 



(12) 



These limits determine the behavior of the polarization 
charge at small and large distances, respectively, which 
we proceed to investigate in more detail. 

The distribution of the polarization charge density in 
real space is given by 



p(r) = Z J ^y(q)n(q)e^ 



(13) 



On the scale of the lattice spacing, a, away from the 
impurity there is a contribution to the screening charge 
that reads 



P(r) 



-Za- <5(r), 



0, 



(14) 



which is valid in the continuum limit (a — ► 0), and comes 
from the linear contribution in (11). In the massless 
situation (m = 0) this localized polarization charge is 
the final result. The finite mass introduces new behav- 
ior, namely additional polarization charge appears dis- 
tributed in space. The full behavior of p(r) is shown in 
Fig. 3 (solid blue line). One clearly identifies two distinct 
regimes whose asymptotic regions are determined by the 
Compton wavelength, the characteristic length scale of 
the problem: 



h 

mv 



1 

m 



(ft 



1) 



(15) 



Q. 





o 


1 E 




o 






O : 




• Exact 


V 




- ocl/r 3 








1 




0.01 



mr 



FIG. 3: (Color online) Polarization charge from Eq. (13) (solid 
line). The effect of additional electron-electron interactions 
(within RPA as discussed in the text), is also shown (dashed 
lines). Inset: magnification of the 1/r 3 behavior at large dis- 
tances. 



Using Eqs. (11,12), for distances much smaller than Ac 
(yet away from the Dirac-delta at r = 0) we find loga- 
rithmic decay, 



p(r) ~ Zam z In I — — 
mr 



ma -C mr <C 1 



(16) 



while at large distances a fast power law emerges 21 

1 



p(v) ~ Zam 



(mrY 



mr 3> 1 



(17) 



These two regimes are evidenced in Fig. 3 by means of 
the log-scale in the main panel [cfr. Eq. (16)] and by 
the fit shown in the inset [cfr. Eq. (17)]. The numerical 
evaluation of Eq. (13) shown in the figure provides the 
crossover between the two asymptotic regimes. 

It is also significant to notice that the two contribu- 
tions — the localized term (14), and the spread tail — 
have different signs: the lattice-scale contribution has a 
screening sign, while the long range tail has a compensat- 
ing, anti- screening, sign. This follows from the fact that, 
per Eq. (12), p(q = 0) = 0, meaning that the total po- 
larization charge Q(oo) = 0, where Q(R) is the vacuum 
charge accumulated within radius R: 



Q{R) - f p(v)dv 

J\t\<R 



(18) 



In fact, the rapid decay of p(r) beyond Ac means that 
most of the additional (positive) charge, compensating 
Eq. (14), is accumulated between the lattice scale (r ~ a) 
and Ac, in such a way that Q(R> 1/m) m 0. This has 
peculiar consequences for screening: the impurity poten- 
tial is best screened at the shortest distances (r ~ a), 
screening weakens between a < r < Ac, and is essen- 
tially absent beyond Ac- Thus the impurity charge re- 
mains unscreened at large distances, as expected for an 
insulator. 
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FIG. 4: (Color online) Exact diagonalization plots of Q(R), 
Eq. (18), in the subcritical regime (we set a — 1 in the plot), 
(a) Finite size gap only, with Ac — 50. (b) Explicit mass 
mv 2 = O.li, leading to Ac — 15. In both cases, the appar- 
ent inflections at R ~ 55 are due to boundary effects of the 
finite system. The top insets compare the maximum of Q(R) 
(Qmax) obtained in the lattice with the value expected from 
Eq. (14). 



IV. INDUCED CHARGE ON A FINITE 
LATTICE 

To confirm the applicability of the above results to the 
real problem on a lattice, and to dismiss possible regular- 
ization issues, we have investigated via exact diagonaliza- 
tion the corresponding tight-binding model for graphene, 
where the lattice appears naturally. The induced charge 
density can be straightforwardly obtained with the aid of 
the exact wave-functions: 



P(r) 



J2 *L(r)^Mr) 



E<- 



E<-i 



^ f (r)^(r) 



(19) 



Summation over the two spin components is also assumed 
(leading to an additional factor of 2). Here 'F^. are the 
wave-functions without external field [Z = 0). Rather 
than address the induced charge itself, it is more con- 
venient to consider Q(R), as defined in Eq. (18). This 
quantity is already averaged over all directions and is 
thus smoother and more appropriate for a direct numer- 
ical comparison. 



We have studied two cases: (a) tight-binding model on 
a finite lattice with 124 x 124 sites, without explicit mass 
term, and (b) the same system with an explicit mass of 
mv 2 = O.li, where t is the hopping parameter. These 
two cases allow us to address the two asymptotic regimes 
in Eqs. (16) and (17). In case (a), although m is ex- 
plicitly zero, there is an effective gap due to the finite 
size of the system, 2mv 2 = A ~ 0.06i, and an effective 
Ac — 50a. 7 Therefore, there is an appreciable region sat- 
isfying a < r < Ac*. The calculated Q(R) for this case 
is shown in Fig. 4(a). Its behavior consists of a sharp 
increase for R ~ a and a subsequent slow decay up to 
Ac- This decay follows the law —Q(R) cx (const. — R 2 ), 
as one expects from Eq. (16). Unfortunately, for our sys- 
tem R = 55a w Ac is the largest distance from the im- 
purity that is free from boundary effects, and one cannot 
comment on the crossover at larger distances. In order 
to address that limit we look at case (b), for which Ac 
is much smaller (Ac ~ 15a); our results are shown in 
Fig. 4(b). In effect, we obtain Q{R) oc in the region 
r ~ ^Ci m accordance with the result in (17). In the 
lower inset of the figure we show the r -3 decay of the ac- 
tual induced charge on the lattice, for a particular value 
of the coupling. The smallness of Ac in this case means 
that the intermediate, logarithmic, regime is inaccessible. 
The analytical behavior thus stands in the lattice prob- 
lem, with qualitative and quantitative agreement, even 
for the case when the gap is due to the finite lattice size. 



ROLE OF ELECTRON-ELECTRON 
INTERACTIONS 



Electron-election interactions can influence the be- 
havior described above. Although questionable on ac- 
count of the strict zero carrier density, one can perform 
re-summation of polarization loops within the random 
phase approximation (RPA) . For weak interactions RPA 
is expected to work quantitatively well, with deviations 
increasing as the interaction a increases. 22 RPA amounts 
to the substitution II(q) -> II(q)/[l - V(q)II(q)] in 
Eq. (13). The outcome of this procedure is given in Fig. 3 
for different values of the interaction, being clear that the 
result is a small decrease in the coefficient of the log in 
Eq. (16). 

In addition, a qualitatively important effect arises from 
self-energy corrections to the polarization. We evaluate 
the self-energy at Hartree-Fock level, which was first done 
for the massless case in Ref. 23. In the massive case we 
obtain a velocity rcnormalization of 



v(q) 




(20) 



where A ~ 1/a is the ultraviolet cut-off, and the above 
expression is valid for q, m <C A. The mass is also rcnor- 
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malized to a larger value in: 















v m)} 


, TO 











<A. (21) 



It is interesting to note that the logarithmic mass renor- 
malization formula in graphene (21) is similar to the well- 
known expression for the electromagnetic mass (account- 
ing for radiative corrections) in 3D relativistic QED. 24,25 
From (20), for A ^> q 2> to, one has v — > v(q) = 
v[l + jln(A/g)], which leads to the "running" of the 
coupling constant: ct(q) = e 2 /v(q). Consequently, ex- 
panding v(q) in powers of a leads, perturbatively, to an 
additional piece in the polarization charge: 

_ , . Zo? 1 , , . 

Sp(r) ~ a < r < 1/to . (22) 

16 r z 

This is the perturbative limit of the effect, first discussed 
in Ref. 3. 

We conclude that electron-electron interactions affect 
somewhat the above discussed behavior at the scale Ac, 
but do not change the picture qualitatively. This can 
be credited to the mentioned absence of screening, and 
the fact that interactions do not generate an additional 
length-scale. Furthermore, it seems that in current ex- 
periments a <C 1, making the interaction corrections 
paramctrically small. 

VI. CONCLUSIONS 

We have found that in the presence of a finite mass gap 
the polarization charge, induced by a Coulomb impurity 
in graphene, has a non-trivial behavior as a function of 
distance. While at zero mass it is concentrated at the 
impurity site, at finite mass it is distributed mostly up 
to Ac* = 1/to, with an additional power law tail beyond 



that. Unlike the massless case, the total vacuum charge is 
zero since the finite mass "pulls" a compensating charge 
from infinity to a distance ~ Ac, and the impurity charge 
remains unscreened beyond this scale. 

In relativistic QED 25 the polarization charge at short 
distances has anti-screening sign (enhances the poten- 
tial), while the compensating charge is distributed up to 
the Compton wavelength of the electron. This behavior 
arises from the running of the charge e 2 (r) in QED. In 
non-relativistic graphene (where the charge is not renor- 
malized), the situation is reversed as shown by the sign 
of Eqs. (14) and (16), (17). 

In the experiment of Ref. 12, where the spectral gap 
is A sw 0.26 eV, we have Ac ~ 30a « 4 nm and the 
behavior we predict in this work should be observable 
if an external ion is present and generates the discussed 
charge re-distribution. In practical terms, in order to 
detect the density variation the chemical potential might 
have to be at or above the value of the gap, > to. Our 
results will then be strictly valid only for \fx\ > to, where 
the screening length (determined by the occupation of 
the conduction band) remains large and well separated 
from the scale Ac ■ Although we strictly have a hyperbolic 
band, when \fi\ > to we can resort to the behavior of a 
parabolic band in 2D. Screening in this case is somewhat 
peculiar 26 due to the finite density of states at the band 
edge, N(fi = to) cx to. Wc expect that this could in 
fact facilitate detection of density variations via STM 13,16 
compared to the massless case in graphene. 
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